Photoacoustic spectrum analysis for spherical target size and optical property determination: A feasibility study

The photoacoustic signal generated by an optically absorbing target is determined by the spatial profile of absorbed optical energy within the target. The analysis of the time profile and frequency content of the signal enables the recovery of the geometry of the object, as well as information about the optical properties. The photoacoustic response of spheres with a homogeneous absorbed optical energy profile is well described, and it is known that the width of the photoacoustic pulse is determined by the diameter of the sphere and its sound speed. In practice, the optical attenuation coefficients within the sphere will result in an inwardly decaying fluence profile leading to a similarly decaying absorbed optical energy profile. Further, the optical attenuation coefficients may be inhomogeneously distributed in the sphere. The implication for both cases is that the existing model for spheres does not fully apply. In this work, we developed analytical expressions for the photoacoustic time traces and amplitude spectra generated by a sphere with absorbed optical energy only in a spherical shell, and by a sphere with an inwardly decaying optical energy profile. Numerical simulations and experiments were conducted on these two imperfect sphere types. Fitting our model to the simulated or measured spectra allowed us to test our model’s ability to extract the sphere size and optical properties. We found that the radii can be recovered with high accuracy, even when the frequency response of the detector recording the photoacoustic pulse is not precisely known. The model was found to be less sensitive in recovering the optical attenuation coefficient, but it is feasible when the detector’s frequency response is well known.


Introduction
Photoacoustic (PA) images are usually created based on tomographic principles, where PA time-traces recorded at several detection points are used to reconstruct images of the optical absorbers. Here a higher density and a larger spatial spread of detection points is known to improve the image quality [1][2][3]. Acquisition of three dimensional PA images with good imaging quality therefore results in high system costs when high numbers of detectors are used, or relatively long measurement times if the detection aperture is synthesized by scanning of a smaller number of detectors.
The PA behaviour of sources with a simple geometry, like spheres, cylinders and slabs, is well known and is described in detail in literature [4,5]. From these, we know that the PA signals contain information about the dimensions of the object. Spheres are, for example, known to generate N-shaped time-traces, with a pulse width that is proportional to the sphere diameter, and with a time of arrival that is proportional to the distance between the source and the detector [5,6]. The amplitude spectrum of such an N-shaped pulse shows an oscillatory behaviour, the PA amplitude, and the reduced scattering coefficient ( ′ ) from a combination of the two. This gave good results in slabs, and should potentially also be applicable in targets with other geometries.
A model describing the PA response of spherical targets with a non-uniform initial pressure profile can potentially find use in various applications. In the biomedical field, at the nano-or micron-scale, one could think of localizing and characterizing exogenous PA or ultrasound (US) contrast or therapeutic agents such as microdroplets [9,13], microparticles [14] or nanostructures [9,13,15] which possess optically absorbing cores or shells. One can also think about macro-scale applications in for example tumour treatment monitoring. Tumour masses can exhibit non-uniform intrinsic optical absorption profiles due to non-uniform vascular distributions. Some tumours types are known exhibit a necrotic core encapsulated by a viable vessel rich rim [16][17][18], while smaller or other tumour types demonstrate less necrosis and have vessel growth throughout the entire tumour [19]. Solid and rim like appearances have also been observed in mice studies with PA contrast agents to monitor drug delivery or treatment progression, [20] where tumours presented high PA signal at the rim shortly after injection, but the intensity in the core increased after 24 h.
In this work, we make the first steps in the development of a model that describes the PA response of spheres with a non-uniform initial pressure distribution. Analytical expressions for the PA time-traces of spherical shells and optically homogeneous spheres with an inward decaying initial pressure profile, and the Fourier transforms thereof, were derived. Amplitude spectra were focused on since time-traces from different objects can be difficult to discriminate due distortions in signal shape as a result of detector frequency band limitations [21], while the locations of maxima and minima in the frequency spectrum, that relate to the sphere size, remain unchanged. This feature together with the large number of peaks in the frequency spectrum are expected to enable a robust and accurate determination of the sphere's dimensions and properties.
The ability of our model to recover the inner and outer radii of the spherical shell, and the radius and the optical attenuation coefficient of the homogeneous beads was investigated with a simulation study in which we tried to fit our models to the spectra from a simulated experiment. To test the model in practice, we measured the PA response of two optically absorbing spherical shells with known radii, and of two absorbing homogeneous spheres with known radii and optical properties. We fitted our model to the measured spectra to recover those parameters. Both the simulations and experiments showed that the model is able to accurately measure the sphere radii. An accurate estimation of the optical attenuation coefficient is more difficult, but simulations show that this is possible when the detector frequency response is known accurately.

Theory
To derive expressions for the PA behaviour of spherical shells and spheres with an inwardly decaying pressure profile, we first have to understand the PA behaviour of a homogeneous sphere. The acoustic transient emitted by a spherical object with a homogeneous initial pressure profile recorded at a distance r from the centre of the sphere has a typical N-shape [5,6,22] as shown in Fig. 1. The width of the time-trace is equal to the time the wave requires to traverse the sphere diameter, the amplitude scales with the absorbed optical energy, and depends on the thermoelastic properties of the object. The amplitude spectrum of an N-shaped pulse shows an oscillatory behaviour, with peak amplitudes decreasing with frequency [8] and with the maximum of the spectrum relating to the sphere size ( ) and sound speed in the sphere ( ) as = 0.33 ∕ [7]. The general expression for the PA transient originating from a spherical target yields [22], Where is the optical absorption coefficient, the fluence, the dimensionless Grüneisen coefficient, the Heaviside step function and the retarded time defined as: In this expression, it is assumed that the sphere is in free suspension in a surrounding medium with the same density and a speed of sound . The presence of impedance mismatch with the surrounding medium would modify the time-of-flight for the acoustic transient to reach the detector and create a set of echoes due to the partial reflection at the sphere boundaries.

Spherical shell
The spherical shell model can be described as a sphere with a transparent core of radius surrounded by an optically absorbing shell of a certain thickness such that the outer radius is . The layer's optical properties are such that the fluence distribution is considered homogeneous within the shell resulting in a uniform initial pressure in the shell. Considering this, we can reasonably assume that the initial pressure field is mathematically equivalent to the difference between two spheres both with homogeneous initial pressure, and radii and ( < ) As illustrated in Fig. 1, this transient consists of a positive and a negative peak separated in time, by the time the PA wave requires to travel through the non-absorbing core. The Fourier transform of the time-dependent expression gives the amplitude spectrum of the PA signal coming from a spherical shell as: with = 2 , = 2 , = arctan ( ) , = 2 √ 2 + 2 2 and = ( − 1). The spectrum is the sum of two periodic terms presenting different modulation frequencies: a ''slow'' modulation in − combined with a ''fast'' modulation in + , and an intermediate in . The spectrum can be predicted to appear as a periodically modulated sine with a global decay in amplitude due to the presence of a pre-factor with inverse dependence on frequency (see also Fig. 1).

Sphere with inward decaying fluence
Here we consider a sphere with radius and homogeneous optical properties such that there is an inwardly decaying fluence distribution resulting in an inwardly decaying initial pressure distribution. Considering a homogeneous illumination on the sphere surface and the incidence angle of photons at the surface random, each point of the surface of the bead acts as an isotropic point source which can be expressed as: [23] Here is the distance from the source, 0 the amplitude of the point source, = 1 3( + ′ ) and = √ . The total fluence ( ) at a given location inside the sphere is the integration over the sphere surface M. Dantuma et al. We know that for instantaneous illumination of an object, the PA pressure signal can be expressed as: [24] (⃗, ) = 1 with the initial pressure at location ⃗′ For a compact notation, the time window is restricted to . The initial pressure locations contributing at the instant to the pressure signal at the detector location ⃗ = (0, 0, ) are forming a spherical cap following 2 + 2 + ( − ) 2 = ( ) 2 inside the emitting sphere. Considering the bead surface defined as 2 + 2 + 2 = 2 , the spherical cap area weighted with the initial pressure values allows to write with ( ) 2 = 2 + ( ) 2 − 2 √ ( ) 2 − ′2 and the maximum radius of the projection of the cap on the xy-plane Finally the acoustic transient can be expressed as The corresponding spectrum, normalized by the prefactor 0 4 , can be expressed aŝ with = 2 ∕ corresponding to the pulse total duration.

Simulation study
Simulations on spheres with different dimensions and were performed to investigate the applicability of our model. With sphere parameters given as input, the generation and propagation of the PA wave was modelled with the k-wave toolbox [25]. Simulations were performed on an isotropic grid with a 40 μm spacing leading to a maximum supported frequency of 18.5 MHz. The sound speed and the density were set to 1000 kg/m 3 and 1482 m/s respectively throughout the entire grid. As a first test, a perfect situation was modelled to verify the analytical expressions for the amplitude spectra of the two sphere types (Eqs. (5) and (13)). The time traces were recorded with an ideal point detector (no bandwidth limitations) that was located at 1.25 cm from the centre of the sphere. Two spheres were considered. The first was a spherical shell with = 4.8 mm and = 4.2 mm, and the second a sphere with = 4.8 mm. The second sphere was imparted with a = 5 cm −1 , which results in an initial pressure profile decaying towards the centre of the sphere following Eq. (7).
In practice, the ultrasound detector behaviour is not ideal, with a frequency dependent and band-limited sensitivity. Therefore, for experiments or modelling studies, one has to assume that the ultrasound detector has a non-ideal frequency response function. This function can be known from data sheets provided by the manufacturer or can be obtained from earlier experiments using known characterization methods. We call this the ''assumed frequency response''. However, the quality of a detector may have degraded in time, a detector may be used with different electronics than the earlier characterization measurements, or errors may have been made during the characterization measurements. The actual frequency response at the time of the physical or digital experiment, we call the ''true detector response'', in contrast to the previously characterized frequency response or the ''assumed frequency response''.
A more realistic set of simulations was performed to investigate the applicability of our model in practice. This true detector response was modelled as a sum of three partially overlapping Gaussians covering a wide band from 0-10 MHz, as in an expected real hydrophone response. Time traces following Eq. (3) and (12) were filtered with this true detector response, and random noise with a maximum amplitude of 0.2 times the signal's standard deviation was added.
To account for the discrepancy between the true and assumed frequency response, we multiplied the analytical spectra with the assumed detector response, and fitted this to the simulated spectra by minimizing the squared error between the two spectra for sphere parameters. The simulated time signals, the true detector responses and the assumed detector responses are shown in Fig. 4. This simulation was performed for 100 different randomly chosen sphere parameters, where always had to be larger than for the spherical shell and varied between 1 and 10 mm. For the sphere with inward decaying fluence, was varied between 0 and 100 cm −1 , and 1/ always had to be larger than /2 to ensure that the pressure decay was not too localized to ensure sufficient frequency content to fall within the bandwidth of the simulated detectors. The true detector response was different for each simulation, with random variations between certain bounds in bandwidth and centre frequencies for the three summed Gaussians. The centre frequencies were varied within 0.5 ± 0.3, 2 ± 1.75 and 5 ± 1 MHz, while their respective bandwidths were varied within 150±30, 100±30 and 80 ± 30% respectively. The imperfect knowledge of the detector response used for the fitting was modelled by adding random errors up to ± 10% to the centre frequencies, amplitude and bandwidths of the Gaussians constituting the real detector response.

Bead fabrication
To verify the applicability of the model in practice, PA measurements were performed on two sets of beads. The first set contained two spherical shell-like beads, each having a different radius and shell thickness. The second set consisted of two homogeneous beads representing spheres with inward decaying fluence, each with a different radius and a different optical absorption coefficient.
The beads were fabricated from sodium alginate (SA) hydrogel (Sigma Aldrich), consisting of the two monomers -L-guluronic acid (39 v/v%) and -D-mannuronic acid (61 v/v%) [26]. A 3 w/v% SA solution was prepared and divided into three batches. Gold nanosphere dispersions (AuNS, 13 nm diameter) were also prepared according to the protocol described in Ref. [27]. The AuNS dispersions were added to two of the batches to end with a 2 w/v% and a 1.7 w/v% SA solution corresponding to concentrations of 2.74 × 10 10 and 4.3 × 10 10 particles/ml. No gold nanospheres were added to the third batch. The SA-AuNS mixtures were optically characterized using transmission measurements in a spectrophotometer (Shimadzu, UV2600, Japan). Beads were produced by the drop-by-drop deposition of the mixtures from 8 mm diameter pipettes into a 0.7 M CaCl 2 cross-linking solution, from a deposition height of approximately 20 cm (see Fig. 2A). The pipette diameter and the deposition height influenced the bead diameter and sphericity respectively. After around 15 minutes the beads sank to the bottom of the CaCl 2 bath, indicating successful gelation.  To create spherical shells, blanco 3 w/v% SA gel beads were first produced and then coated with the absorbing SA-AuNS gel in a subsequent step. The cross-linked blanco beads were immersed in the SA-AuNS solution, assuring total coverage of the bead. Re-deposition of these beads in the CaCl 2 bath gelated the SA-AuNS layer around the SA bead. Different shell thicknesses were obtained by varying the immersion time of the bead in the SA-AuNS solution. Fig. 2A illustrates the preparation process.
Several beads were produced of which four beads were selected based on their radii. To centralize a bead inside the measurement setup, each selected bead was embedded into an agarose ball (see Fig. 2B for a schematic of the protocol). The agarose ball was made by casting an aqueous solution of agarose (3 w/v%), into a spherical mould, which consisted of two 3D printed hemispherical moulds (5 cm diameter) that could be combined. One of the hemispherical moulds was filled with the agarose solution and a bead was positioned on the surface when the agarose had partially solidified. After full solidification, it was flipped and placed on top of the other half sphere which had been filled with the agarose solution, allowing the two halves to stick together. The mould was removed after the solidification of the agarose. After the PA measurements, the SA beads were removed from the agarose and cut into two halves to measure their radii with a microscope.

Experimental setup and measurement procedure
The agarose balls carrying their respective gel beads were positioned centrally in a 3D-printed holder (shown in Fig. 2.E&F). The PA signal generation was achieved by coupling 10 ns 532 nm laser pulses generated by a Quanta-Ray Pro-250 Nd:YAG laser (Spectra Physics, USA) with a 10 Hz repetition rate into a multi-output fibre bundle (CeramOptec GmbH, Germany). The pulse energy was distributed equally over the outputs with an average pulse energy per fibre bundle output of 15 mJ. Considering the 10 ns pulse width and 4 mm fibre diameter, one can calculate that the power density at each fibre output yielded 1.5 MW cm −2 . The fibre bundles were placed in the holder in direct contact with the agarose ball from 6 sides to aim for homogeneous illumination of the bead. A PVDF hydrophone detector (Precision Acoustics, United Kingdom) with a bandwidth from 0.1-20 MHz was used to detect the acoustic transients from the beads, instead of the typical bandlimited detector used in PA imaging. The reason for this choice was to allow fitting of the amplitude spectra over the broad range since this is a feasibility study of the proposed approach. The active element (1 mm diameter) of the hydrophone was positioned a few mm from the agarose ball surface. The hydrophone was connected to a submersible pre-amplifier (Precision Acoustics, UK) and a DC coupler (Precision Acoustics, UK) followed by a desktop PC with built-in digitizer (DP105, Acqiris, Switzerland) at 500 Ms/s to save the time traces. Averaging was used to increase the signal-to-noise ratio of the recorded signal.

Simulation study
The unrealistic simulations on two beads using an ideal ultrasound detector have demonstrated that the analytic expressions for the amplitude spectra are correct. Fig. 3 presents the initial pressure distribution slices of the simulated spheres, the recorded time traces and corresponding amplitude spectra. The time traces and amplitude spectra are overlaid with the plots from the analytical expressions (Eqs. (3), (5), (12) and (13)) using the known sphere dimensions and/or attenuation coefficient as input. Despite the noise on the time traces, likely originating from the finite grid size in the simulation, a smooth amplitude spectrum is obtained for both cases and they are in good agreement with our analytic spectra. Fig. 4 shows the simulated time traces and frequency spectra for two out of the 200 simulated experiments. The fitted spectra together with the found radii and/or attenuation coefficient, and the used detector responses are also presented in the Figure. Both for the spherical shell as for the decaying fluence, our model was able to recover the radii very accurately. For both spheres, was found to be equal to exactly the set size.
was slightly overestimated by 0.002% and was overestimated by 13%.
Linear regression analysis was performed on the found radii and for the 100 spherical shells and 100 spheres with decaying initial pressure profiles (see Fig. 5). R-squared values of 99% for the sphere radii were found while withstanding errors in the knowledge of the detector responses. The model seemed to have more difficulties in recovering , especially for higher . An R-squared value of 76% was found. Simulations with a known detector response lifted the Rsquared value significantly to 98%. Experiments with a well-calibrated detector could therefore be accurate in recovering .   Fig. 6 shows photographs of the selected beads and the absorption spectra of the two SA-AuNS dispersions used in making the beads. The left panel of Fig. 7 shows the time traces of the measurements on the beads. Although the SNR of the signals is limited, two peaks corresponding to the inward and outward propagating PA waves can be observed in all four signals. The signals show some deviations from the theoretical models, (compare with Fig. 4) as a result of band-limited ultrasound detection, internal reflections and imperfect sphericity of the beads. The model did converge for both spherical shells meaning that fitting radii were found. The middle panel of Fig. 7 shows the squared errors between the experimental spectrum and the theoretical spectra for varying radii. The cross cursor highlights the found minimum in each case. The spectra corresponding to the parameters of the found minima are plotted in the right panel of Fig. 7 together with the experimentally measured spectra. For the first 2 MHz of beads 1 and 2, the peak locations in the fitted amplitude spectra match well with the locations of the measured spectrum. For the higher frequencies, the fit becomes less accurate. The radii estimated with our model, are = 2.66 mm and = 2.28 mm for bead 1, and = 3.33 mm and = 2.40 mm for bead 2. For both beads, the estimated radii match very well with the sizes observed in the microscopy images which yield 2.63 mm ( ) and 2.07 mm ( ) for bead 1 and 3.26 mm ( ) and 2.43 mm ( ) for bead 2.

Experimental study
Fitting our model to the measurements on the homogeneous spheres (beads 3 and 4) was more complicated, as was already realized from the simulation study. The experimental decays can be observed in the two peaks in the time traces which also results in exponentially decaying amplitude spectra. The model was again able to accurately recover the bead radii. Radii of 2.55 mm and 2.02 mm were found for bead 3 and bead 4 respectively, which were measured 2.52 mm and 2.10 mm in the microscopy images. The estimated attenuation coefficients were however highly overestimated. Coefficients of 19.1 cm −1 and 47.4 cm −1 were found, compared to the ground-truth values of 3.73 cm −1 and 5.78 cm −1 . The vertical smeared-out lines in the squared error figures (middle panel) confirm that the sensitivity of the model in estimating is poor compared to the accuracy in radius estimation.

Discussion
The model as presented in this chapter has a few limitations that should be kept in mind. Firstly, the bandwidth of the generated PA signals should fall within the bandwidth of the detector. When using band-limited detectors, the frequency coverage must be enough to allow identification of either the double modulation in the shell case or single modulation and decay for the homogeneous sphere case. For example, a low-frequency detector used to measure PA signals emitted by a shell structure with a thickness smaller than the main sensitivity wavelength of the detector does not resolve the time trace enough to visualize both modulations in the spectrum. A higher frequency detector with broader absolute coverage may be more appropriate if nothing is known a priori about the target.
As also seen in this work, the calibration of detectors needs to be accurate, as the frequency response is implemented into the fitting. In particular, for the homogeneous sphere case, where the slope of the envelope of the spectrum informs on the optical properties, accurate frequency sensitivity estimation is essential.
The radii and can only be estimated accurately if the peaks and valleys in the spectra are resolved well. Zero padding of the time traces can suffice, with the consequence of increased computation times. This is usually not a problem for a single target study. But more intensive applications with iterative computation for multiple targets of the corresponding model and detector response can become time-consuming. The signal-to-noise ratio also impacts the quality of the obtained spectra. It can be experimentally improved with simple implementations such as averaging or more efficient delivery of the laser light to the bead or more sensitive detectors.
Then there are also some limitations of the geometrical parameters of the spherical shell. If the considered shell thickness is large compared to the inner core, the proposed approach is prone to failure. In the analytical formulation, the argument of the product of sine and cosine tends to the same value. The double periodicity becomes subtle and in addition to accurate detector calibration, a finer frequency step size is needed to obtain an accurate fit.
In the current definition of the shell model we suppose homogeneous initial pressure in the absorbing layer. This hypothesis holds when the shell is 'optically-thin' where 1∕ ≫ − , meaning that most of the light is transmitted by the shell and contributes to the fluence distribution on the other side of the shell. If 1∕ ≫ − , the fluence will decay within the shell and the model for homogeneous beads with an inward decaying fluence might be more appropriate to use. Thermal and stress confinement is assumed in our model but is only valid for macro and micro-scaled absorbers excited by nanosecond light pulses. For nano-scale absorbers, the thermal and stress confinement can be violated, meaning that the model loses accuracy. It is known from experimental and theoretical work that PA signals from gold nanoparticles, due to non-conformity of thermal confinement, are caused by heat transfer across the nanoparticle-medium interface and subsequent thermal expansion in the surrounding medium. The temperature dependence of water thermal expansion can provide amplification of the PA signals [28][29][30]. According to recent theoretical work [31], the non-linear effects in the proximity of the plasmon resonance peak, can increase the PA signal by more than 50%. This can be an additional contribution to the overestimation in observed in our case, but cannot be confirmed in the scope of the present work.
The model can be made more accurate by incorporating more physical properties of the media. Currently, the model assumes an acoustically homogeneous medium and does not account for impedance mismatches that give rise to acoustic reflections between e.g. the sphere and the surrounding medium or the shell and the core of the sphere. In our experiments, hydro-gels in water were used to have minimal impedance mismatches and thereby minimal reflections. In possible future applications, reflections will likely take place which may make the model fitting more challenging. Also, the frequency-dependent acoustic attenuation inside the bead and in the water are not included in the model. High frequencies attenuate more than low frequencies and can thereby complicate the fitting of the . Implementation of more physical properties, like the ones mentioned above, are expected to increase the model's ability in property determination in experimental situations, but will also make the model more complex. The level of complexity that is required differs per application and is something to be researched in the future.
Experimental imperfections like irregular distributions of optical properties within the shell or sphere, inhomogeneous illumination or imperfect sphericity of the target, could complicate the fitting of our model. Although the SA beads used in our experiment also had imperfect sphericity the model was still able to accurately estimate the radii. This provides the first indications that the model can be robust against imperfections. Further research with simulations on imperfect spheres should point out how these imperfections exactly influence the outcomes and for what imperfections the model is still able to find accurate solutions. Alternatively, some of these imperfections, like the imperfect sphericity could potentially also be implemented in the model in the future.
The presented approach and the feasibility study is meant specifically for measuring the size/optical properties of spherical shaped structures that have been targeted by PA or US contrast agents. These contrast agents extravasate and accumulate either passively or actively at a tumour site. However, the distribution pattern and leakiness of blood vessels in tumours is extremely heterogeneous. Specifically, the microvascular density is high at the invasive edge of the tumour, but sometimes the tumour centre is malperfused, hindering delivery of the contrast agent to the interior region. This is the reason for studying a rim-distribution of absorption (spherical shell), in addition to the homogeneous distribution with realistic radial-decay of fluence.
The method introduced in this work can form the basis of translation to cylindrically-shaped objects. New analytical expressions will have to be derived based on the PA behaviour of cylindrical absorbers. Similar methodologies as proposed can be used, in which case the estimation can also be attempted for cylindrical or line sources and then further on for blood vessels. For those cases, the PA spectrum will be dependent on the orientation of the absorber with respect to the US detector [32]. A caveat is that in vivo PA signals will be the integral of multiple sources within the acceptance angle of the detector, so that extraction of the frequency components from specific structures and locations will be challenging.

Conclusion
In this work, we have proposed an approach to estimate the radii and optical attenuation of spherical objects from a single-point measurement of the PA signals produced. We discuss two cases -that of a sphere with pressure distribution in a shell due to optical absorption at the surface region only, and that of a sphere with inwardly decaying pressure distribution following the fluence profile in an homogeneously absorbing object. We present analytical expressions describing the time traces and amplitude spectra of the PA response of the two cases. The ability of this model in extracting the radii of spherical shells, and radii and optical attenuation coefficients of spheres was first investigated with a numerical experiment where simulations were performed to mimic an experiment. This study showed that the model was able to recover the radii without problems, but that accurate knowledge of the detector frequency response is needed to also achieve accurate estimations of the optical attenuation coefficient. As a next step, an experiment was performed on two spherical shells and two homogeneous spheres made from sodium alginate doped with gold nanoparticles. Experiments on those beads showed comparable outcomes to the simulation study. The radii were recovered accurately, but the estimation of the optical attenuation for the homogeneous spheres was inaccurate. In future studies, a broadband detector with a more accurately known detector response should be used to further investigate the ability of the model to assess the optical attenuation coefficient.

Funding
This work was part of the European Horizon 2020 PAMMOTH project under grant agreement No 732411, an initiative of the Photonics Public Private Partnership.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.
Maura Dantuma (Ph.D.) earned her Bachelor degree in the biomedical photonic imaging group at the University of Twente, The Netherlands, on the subject of contrast agents for photoacoustic imaging in tissue engineering. During her studies she developed interest in biomedical optics and acoustics. She finished her master degrees in Applied Physics and Biomedical engineering in the Physics of Fluids group, University of Twente, where she investigated acoustic droplet vaporization of free space perfluoropentane droplets. For her Ph.D. degree she was part of the Horizon 2020 PAMMOTH project, where she focused on the development, optimization and testing of a multi-wavelength photoacoustic-ultrasound breast tomography system. After a successful completion of the Ph.D. in the Multi-Modality Medical Imaging (M3I) group at the University of Twente, she started as a post-doc in the same group, where she continued working on photoacoustics.